TD2 - 20/04/21 - Simulation et Monte Carlo

• Exercice 4 : Loi géométrique
• Exercice 5 : Control variates, variables antithétiques, QMC
• Exercice 6 : Importance Sampling
• (Exercice 7 : MCMC)

Exercice 4 : Loi géométrique

• Proposer un algorithme pour simuler une loi géométrique
• Le mettre en oeuvre
• Proposer une méthode pour vérifier les résultats

Méthode 1 :

soit $X\sim\mathcal{G}(p)$ qui suit une loi géométrique de paramètre 𝑝.

$$ \forall k \in \mathbb{N}^*, \ \ \ P(X = k) = (p-1)^{k-1}\cdot p = q^{k-1}\cdot p $$

C'est la loi d'attente d'un premier succès dans une suite d'épreuves répétées indépendantes ayant chacun la probabilité de succès 𝑝

Par conséquent si : $$ (U_i)_{i\geq 1} \sim^{i.i.d} \mathcal{U}(0,1) $$

Alors la variable : $$ Y := \min\{k \in \mathbb{N}^* : U_k < p\} \sim \mathcal{G}(p) $$

Méthode 2 :

code :

proposition de générateur de la loi géométrique

import numpy as np

# méthode 1
def rand_geom(n, p):
    """
    n : taille de l'échantillon
    p : paramètre de la loi géométrique
    """
    rg = np.zeros(n)
    for i in range(n):
        ru = np.random.rand()
        k = 1
        while ru > p:
            k += 1
            ru = np.random.rand()
        rg[i] = k
    return rg

# méthode 2
def rand_geom2(n, p):
    """
    n : taille de l'échantillon
    p : paramètre de la loi géométrique
    """
    ru = np.random.rand(n)
    X = np.log(ru)/np.log(1-p)
    Y = np.floor(X)+1
    return(Y)

def proba_geom(k, p):
    return((1-p)**(k-1) * p)

générons un échantillon de taille (exemple) $n = 10000$ avec le paramètre $p = 0.3$

n = 10000
p = 0.3
sample_geom = rand_geom(n, p)
sample_geom2 = rand_geom2(n, p)

affichons son histogramme ainsi que la vraie distribution de la loi géométrique

import matplotlib.pyplot as plt

plt.hist(sample_geom, color = "steelblue", 
         density = True, 
         bins = int(sample_geom.max()), 
         alpha = 0.5, 
         label = "Notre échantillon (1)")
plt.hist(sample_geom2, color = "red", 
         density = True, 
         bins = int(sample_geom2.max()),
         alpha = 0.5, 
         label = "Notre échantillon (2)")
x = np.arange(1, max(sample_geom.max(), sample_geom2.max()))
y = proba_geom(x, p)
plt.plot(x, y, color = "orange", label = "Vrai proba", marker = "o")
plt.xlabel("k")
plt.ylabel("probabilité")
plt.title(f"histogramme et pdf de la loi géométrique : n = {n} et p = {p}")
plt.legend()
plt.show()

Regardons les deux extrêmes : $p = 0.1$ et $p=0.9$

• $p=0.1$

n = 10000
p = 0.1
sample_geom = rand_geom(n, p)
sample_geom2 = rand_geom2(n, p)

plt.hist(sample_geom, color = "steelblue", 
         density = True, 
         bins = int(sample_geom.max()), 
         alpha = 0.5, 
         label = "Notre échantillon (1)")
plt.hist(sample_geom2, color = "red", 
         density = True, 
         bins = int(sample_geom2.max()),
         alpha = 0.5, 
         label = "Notre échantillon (2)")
x = np.arange(1, max(sample_geom.max(), sample_geom2.max()))
y = proba_geom(x, p)
plt.plot(x, y, color = "orange", label = "Vrai proba", marker = "o")
plt.xlabel("k")
plt.ylabel("probabilité")
plt.title(f"histogramme et pdf de la loi géométrique : n = {n} et p = {p}")
plt.legend()
plt.show()

• $p = 0.9$

n = 10000
p = 0.9
sample_geom = rand_geom(n, p)
sample_geom2 = rand_geom2(n, p)

plt.hist(sample_geom, color = "steelblue", 
         density = True, 
         bins = int(sample_geom.max()), 
         alpha = 0.5, 
         label = "Notre échantillon (1)")
plt.hist(sample_geom2, color = "red", 
         density = True, 
         bins = int(sample_geom2.max()),
         alpha = 0.5, 
         label = "Notre échantillon (2)")
x = np.arange(1, max(sample_geom.max(), sample_geom2.max()))
y = proba_geom(x, p)
plt.plot(x, y, color = "orange", label = "Vrai proba", marker = "o")
plt.xlabel("k")
plt.ylabel("probabilité")
plt.title(f"histogramme et pdf de la loi géométrique : n = {n} et p = {p}")
plt.legend()
plt.show()

Exercice 5 : Control variates, variables antithétiques, QMC

5.1.

code

Simulons notre vecteur gaussien de dimension $d\geq2$. Nous pouvons utiliser l'exercice 3 pour les simuler mais il est aussi possible de les simuler via grâce à la librairie numy : numpy.random.multivariate_normal

import numpy as np

d = 3
n = 5000
mu = np.zeros(d)
sigma = np.eye(d)

rn = np.random.multivariate_normal(mu, sigma, n)

Jettons un coup d'oeil à l'échantillon généré (vous pouvez utiliser le module mplot3d de matplotlib)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(rn[:,0], rn[:,1], rn[:,2])
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$x_3$')
plt.show()

Programmons la méthode Monte Carlo :

def phi(x, c):
    return (abs(np.prod(x)) < c)*1

def monte_carlo_estimate(n, d, c):
    """
    n : taille de l'échantillon
    d : dimension
    c : le paramètre pour la région A
    """
    rn = np.random.multivariate_normal(np.zeros(d), np.eye(d), n)
    I_hat = 0
    for i in range(n):
        I_hat = I_hat + phi(rn[i,:], c)
    I_hat = I_hat / n
    return I_hat

def monte_carlo_estimate_list(n, d, c):
    """
    n : taille de l'échantillon
    d : dimension
    c : le paramètre pour la région A
    """
    rn = np.random.multivariate_normal(np.zeros(d), np.eye(d), n)
    rv = [phi(rn[i,:], 0.9) for i in range(n)]
    I_hat = np.cumsum(rv)
    I_hat = I_hat/np.arange(1,len(rv)+1)
    return I_hat

monte_carlo_estimate(n = 100, d = 3, c = 1)

0.82

monte_carlo_estimate(n = 100, d = 3, c = 0.5)

0.73

monte_carlo_estimate(n = 100, d = 3, c = 0.0001) 
# lorsque c est petit, l'événement devient rare et l'estimation donne presque tout le temps 0

0.0

monte_carlo_estimate(n = 10000, d = 3, c = 0.0001)
# il est possible d'avoir une meilleure estimation en augmentant la taille de l'échantillon

0.0029

• regardons l'évolution de l'estimateur en fonction de $n$ pour $c=0.5$

N = 200

x = np.arange(1, N+1)
y = []
for n in x:
    y.append( monte_carlo_estimate(n, d = 3, c = 0.5) )
plt.scatter(x,y)
plt.xlabel("n (taille de l'échantillon)")
plt.ylabel("$\hat I_n$")
plt.show()

• regardons l'évolution de l'estimateur en fonction de $n$ pour $c=0.0001$

x = np.arange(1, N+1)
y = []
for n in x:
    y.append( monte_carlo_estimate(n, d = 3, c = 0.0001) )
plt.scatter(x,y)
plt.xlabel("n (taille de l'échantillon)")
plt.ylabel("$\hat I_n$")
plt.show()

Nous nous situons dans une situation où l'événement est rare et la méthode monte carlo "classique" n'est pas stable pour $n$ petit.

5.2.

Idée : Prendre la statistique $h(X) = \sum\limits^{d}_{k=1}X^k$ comme variable de contrôle, pour deux raisons :

• $h(X)$ est corrélé à $\{X \in A\}$

• il est facile de calculer l'espérance de $h(X)$ : $\mu = \mathbb{E}[h(X)] = d\cdot \mathbb{E}[X^1] = d\cdot \sqrt{\frac{2}{\pi}}$

Jettons un coup d'oeil sur la forme de la fonction $\varphi$ :

def h(x):
    return(sum(abs(x)))

def control_variates_estimate(n, d, c):
    """
    n : taille de l'échantillon
    d : dimension
    c le paramètre pour la région A
    """
    rn = np.random.multivariate_normal(np.zeros(d), np.eye(d), n)

    fX = [phi(rn[i,], c) for i in range(n)]
    fXbar = np.mean(fX)
    
    hX = [h(rn[i,]) for i in range(n)]
    mean_hX = d*np.sqrt(2/np.pi)
    var_hX = np.var(hX)
    
    # "[0][1]" car la fonction "cov" de numpy donne une matrice et l'élement que nous chercons se trouve dans la coord [0][1]
    cov_hX_fX = np.cov(fX, hX)[0][1] 
    beta = - cov_hX_fX / var_hX
    
    Z = fXbar + beta * np.mean(hX - mean_hX)
    return(Z)

n = 200
d = 3
c = 0.01
Imc = []
Icv = []
for i in range(n):
    Imc.append( monte_carlo_estimate(n, d, c) )
    Icv.append( control_variates_estimate(n, d, c) )

plt.boxplot((Imc, Icv), labels = ["Monte Carlo", "Control Variates"])
plt.show()

5.3.

Idée : Nous pouvons générer

• $X_i = G(U_i)$ par Box-Muller avec les variables uniformes $U_1$ et $U_2$ et
• $-X_i = G(1- U_i)$ avec les variables uniformes $1-U_1$ et $1-U_2$

Exercice 6 : Importance Sampling

6.0.

Programmons tout d'abord la méthode Monte Carlo qui va nous servir de baseline

# monte carlo standard

def phi(x):
    # x : un vecteur de dimension 2
    return np.log(1 + abs(x[0] - x[1]))

def monte_carlo_estimate(n, d = 2, c = 0.5):
    """
    n : taille de l'échantillon
    d : dimension (2 par default)
    c le paramètre pour la région A (0.5 par default)
    """
    
    rn = np.random.multivariate_normal(np.zeros(d), np.eye(d), n)
    index = (abs(np.prod(rn, axis=1)) < c)
    rn = rn[index, :]
    m = len(rn)
    
    I_hat = 0
    if m>0:
        for i in range(m):
            I_hat = I_hat + phi(rn[i,])
        I_hat = I_hat / m
    
    return(I_hat, m/n)

monte_carlo_estimate(100, d=2, c=0.01)

(0.17701302520707451, 0.06)

6.1.

Proposition de simulation :

1. Simuler les $X_i \sim \mathcal{N}_d(0, I_d)$

2. Calculer la moyenne $\mu$ et la matrice covariance empirique $\hat \Sigma$ des $X_i$ appartenant à $A$

3. Faire de l'IS avec la loi auxiliaire $\mathcal{N}_d(\hat\mu, \hat \Sigma)$

Calculer $\hat\mu$ et $\hat \Sigma$ est une mauvaise idée car ces termes sont encore des estimateurs de MC standard sur des événements rares lorsque $c$ devient petit ...

Code :

from scipy import stats

def IS_gauss_estimate(n, d = 2, c = 0.5):
    # etape 1
    X = stats.multivariate_normal.rvs(np.zeros(d), np.eye(d), n)
    index = (abs(np.prod(X, axis=1)) < c)
    X = X[index, :]
    m = len(X)

    # etape 2
    mu = np.mean(X, axis=0)
    sig = np.cov(X.T)

    # etape 3
    Y = stats.multivariate_normal.rvs(mu, sig, n)
    g = stats.multivariate_normal.pdf(Y, mu, sig)
    f = stats.multivariate_normal.pdf(Y, np.zeros(d), np.eye(d))
    index = (abs(np.prod(Y, axis=1)) > c)
    f[index] = 0
    m = sum(index)

    W = f / g
    phi_g = np.array(list(map(phi, Y)))
    I_tilde = sum(phi_g * W) / sum(W)
    
    return(I_tilde, m/n)


print(monte_carlo_estimate(10000, d=2, c=0.01))
print(IS_gauss_estimate(10000, d=2, c=0.01))

(0.30832733369902626, 0.0375)
(0.27251511485356356, 0.8557)

6.2.

On change la loi auxiliaire / de proposition $g$ par la t-distribution qui est plus épaisse aux bords qu'une gaussienne

def IS_t_estimate(n, d=2, c=0.1, df=5):

    rnd = np.random.RandomState()
    # On genere une t de taille (n, d)
    Xg = rnd.standard_t(df, size=(n, d))
    # On calcule la densite de la loi t non normalisee
    g = (1 + (Xg ** 2).sum(axis=1) / df) ** (- (2 + df) / 2)
    # et celle de f (gaussienne tronquee)
    f = np.exp(- 0.5 * (Xg ** 2).sum(axis=1))
    
    index = (abs(np.prod(Xg, axis=1)) > c)
    f[index] = 0.
    weights = f / g
    weights = np.nan_to_num(weights)
    m = sum(index)
    
    phi_g = np.array(list(map(phi, Xg)))
    phibar = (phi_g * weights).sum() / weights.sum()
    effective_n = weights.sum() ** 2 / (weights ** 2).sum()
    
    return np.nan_to_num(phibar), m/n

n, d, c = 10000, 2, 0.001

print(monte_carlo_estimate(n, d, c))
print(IS_gauss_estimate(n, d, c))
print(IS_t_estimate(n, d, c))

(0.22413471707847274, 0.0075)
(0.2490232919845994, 0.962)
(0.27342166394769074, 0.9961)

n, d, c = 10000, 2, 0.0001

print(monte_carlo_estimate(n, d, c))
print(IS_gauss_estimate(n, d, c))
print(IS_t_estimate(n, d, c))

(0.1239590618205508, 0.0009)
(0.08881983328552698, 0.9732)
(0.2843876739789567, 0.9996)

n, d, c = 1000, 2, 0.01
N = 50

data = np.zeros((N, 3))
data[:, 0] = np.array([monte_carlo_estimate(n, d, c)[0] for _ in range(N)])
data[:, 1] = np.array([IS_gauss_estimate(n, d, c)[0] for _ in range(N)])
data[:, 2] = np.array([IS_t_estimate(n, d, c)[0] for _ in range(N)])

plt.figure()
plt.boxplot(data, labels=["MC classique", "IS gaussian", "IS t"])
plt.show()

Exercice 7 : MCMC (correction partielle)

7.1.

• Quelle loi reconnaissez-vous ?

7.2.

• Calculer les lois conditionnelles de chaque composante

• Mettre en oeuvre le Gibbs sampler correspondant

code : SOON ....